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ABSTRACT 

In this paper, we study the heating of the magnetized solar chromosphere in- 
duced by the large fraction of neutral atoms present in this layer. The presence of 
neutrals, together with the decrease with height of the collisional coupling, leads 
to deviations from the classical MHD behavior of the chromospheric plasma. A 
relative net motion appears between the neutral and ionized components, usually 
referred to as ambipolar diffusion. The dissipation of currents in the chromo- 
sphere is enhanced orders of magnitude due to the action of ambipolar diffusion, 
as compared to the standard ohmic diffusion. We propose that a significant 
amount of magnetic energy can be released to the chromosphere just by existing 
force-free 10-40 G magnetic fields there. As a consequence, we conclude that 
ambipolar diffusion is an important process that should be included in chromo- 
spheric heating models, as it has the potential to rapidly heat the chromosphere. 
We perform analytical estimations and numerical simulations to prove this idea. 

Subject headings: Sun: chromosphere - Sun: magnetic field - Sun: numerical 
simulations 

1. Introduction 

The lower solar atmosphere — photosphere and chromosphere — is usually treated using 
the approximation of magnetohydrodynamics (MHD). Different dynamical processes, like, 
e.g., wave propagation, formation of magnetic structures, flux emergence, etc., have been 
quite successfully described with this formalism. Nevertheless, because of the rather cool 
temperatures of the photosphere and low chromosphere, the degree of ionization there is very 
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small, reaching values as low as 10~ 4 in the tem perature minimum, a nd remaining always 
well below unity even at larger heights (see, e.g. IVernazza et al.lll98ll ). This fact, together 
with the decrease of collisional coupling with height, leads to a break of the assumptions 
underlying MHD. In the upper photosphere and chromosphere, collisions are unable to couple 
completely the neutral and ionized components. Due to that, several new effects appear, 
such as non-ideal Hall effect and ambipolar diffusion. 

In plasma physics, the term "ambipolar diffusion" refers to the diffusion of positive and 
negative particles at the same rate due to their Coulomb interaction, which maintains the 
charge neutrality at scales larger than the Debye length. In astrophysics, however, ambipolar 
diffusion usually refers to the decoupling of neutral and charged components. Ambipolar 
diffusion causes the magnetic field to diffuse though neutral gas due to collisions between 
neutrals and charged particles, the latter being frozen-in into the magnetic field. The use of 
the same terminology applied to these two different phenomena may lead to confusion. In 
this paper we follow the astrophysical definition of the ambipolar diffusion^! As for the Hall 
effect, it appears as a result of the different drift velocities of electrons and ions, which are 
not equally affected by the presence of neutrals. The spatiotemporal scale over which the 
Hall effect oper ates in a partially ionize d plasma is very different from the case of a fully 
ionized plasma ( IPandey fc Wardldl2008l ). 



Recently, there has been an increasing number of works in the literature showing the 
importance of the deviations from MHD in different situations. The presence of neutral 
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tion between the neutral and charged species causes an increase of the collisional damping of 
MHD waves in the photosphere, chromosphere and prominence plasmas. This can be an im- 
portant way to le ak energy from Alfven and fast m ode waves, causing the dam ping of coronal 



loop oscillations ( IKhodachenko et al.ll2004l 120061 ) and prominence oscillations (IForteza et al. 



20071 ). Under certain conditions, the leakage of Alfven wave ene rgy can become the source 
of instabilities in the medium due to the non-ideal Hall effect (IPandey et al.l 120081 ). The 
excitation rate of Alfven waves by foot-point motions of magnetic structure s significantly 
decreases if the neutral component is taken into account (jVranjes et al.l 120081 ). Another ef- 



fec t is the poss i ble ap pearance of cut -off frequencies for Alfve n and kink waves, as mentioned 



by lSoler et al.l (120091 ) (see, however, iZaqarashvili et al.l 120111 ). 



1 Although the term "ambipolar diffusion" is widely used in the literature, we consider more appropriate 
the term "diffusion by neutrals" , or "neutral diffusion" . In the rest of the paper we will use both terminologies 
interchangeably. 



-3 - 



The plasma par tial ionization i s also important for magnetic reconnection, for the same 
reason as for waves. IZweibell ( 119891 ) has shown that the reconnecti on rate depends strongly 
on the collis ional coupling between ionized and neutral species and iBrandenburg fc Zweibel 
( 119941 Il995l ) concluded that, due to the action of ambipolar diffusion, oppositely oriented 
magnetic field lines can be br o ught sufficiently c lose t o facili tate the reconn e ction. In a series 
of papers, ISakai et al.l (120061 ) . ISmith fc Sakail ( 120081 ). and ISakai fc Smith! (120091 ) applied a 
two-fluid approximation to study the reconnection between two current loops. The two-fluid 
approach allowed these authors to find different temperatures of the ionized and neutral 
species, and give rise to proton heating and jet-like phenomena that, possibly, can be the 
reason for footpoint heating of coronal loops, explosive events and jets in the transition 
region and sunspot penumbrae. 

All these non-ideal plasma effects can also modify the equilibrium balance of photo- 
spheric and chromospheric structures. An interesting mechanism to concentrate kilo-Gauss 
magnetic flux tubes at photospher i c leve l, due to the action of the Hall current, was pro- 
posed by iKhodachenko fc Zaitsev! (120021 ) . At the chromosphere, the order of magnitude 
stronger dissipation of currents perpendicular to the magnetic field, compared to that of 
longit udinal currents, wa s found to facilitate the creation of potential force-free field struc- 
tures (lArber et al.l 120091 ) . In prominences, the frictional force generated due to collisions 
between neutrals a nd ions was found to play a role in supporting their structure against the 
gravitational force ( {Gilbert et al.l 120021 ). 



Another phenomenon affected by n on-ideal plasma effects is magnetic flux emergence 
( iLeake &: Arberll2006l ; lArber et al.ll2007l ). These authors showed that the amount of emerged 
flux is greatly increased by the presence of a diffusive layer of partially ionized plasma. 
Additionally, including neutrals removes the non-realistic lifting of low-temperature plasma 
into the lower corona and also avoids the subsequent Rayleigh- Taylor instability. 

All the examples given above indicate that incorporating non-ideal plasma effects is 
essential and can lead to new insights into the physics of many phenomena taking place in 
the solar photosphere and chromosphere. Despite this increasing evidence, we are far from a 
complete understanding of the influence of these effects. The aim of the present paper is to 
investigate the consequences of the Joule dissipation of electric currents into the heating of the 
magnetized chromosphere. In the presence of neutrals, the ambipolar (or neutral) diffusion 
is expected to be orders of magnitude large r than the classical ohm i c diffusion, leadin g 
to Reynolds number around unity (see, e.g., IKhodachenko et al.ll2004j ; lArber et al.l 120071 ). 
Diffusion by neutrals creates, in addition, anisotropy in the proper ties of the plasma since 
currents perpendicular to the magnetic field dissipate much quicker ( lArber et al.ll2009l ). Due 
to that, one might expect important plasma heating. This work was inspired by the recent 
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research by |Pe Pontieu fc Haerendell (119981 ) , IJudgd (120081 ) and iKrasnoselskikh et al.l ( I2010l ) , 
pointing in the direction that current dissipation, enhanced by the presence of neutrals in a 
plasma not entirely coupled by collisions, can play an important role in the energy balance 
of the chromosphere and above. In the rest of the paper, we will show analytical calculations 
and numerical simulations in different models (from purely academic to more realistic) of 
the amount of heating of chromospheric flux tubes due to the ambipolar diffusion effect. 



2. Equations for partially ionized plasma 

Depending on the temporal and spatial scales under study, several approaches can be 
adopted to describe the multi-component solar plasma. Without going into a kinetic de- 
scription (impossible to use in the photosphere and lower chromosphere because they are 
too dense), a detailed approach would consist in describing the three components — neutrals, 
ions and electrons — by a system of separate coupled equations (three-fluid approximation). 
However, when the collisional coupling is large, a less detailed two-fluid (adding up equations 
for electrons and ions) or single-fluid (adding up all three equations) approximations can be 
used. Here we adopt the latter approach. The derivation of the conservation equations 
for in dividual species in this case can be found in, e.g. 



Braeinskii 


(1965 


) and 
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(119861 ). For simplicity, we consider a hydrogen plasma with three components: hydrogen 
ions (i), neutral hydrogen (n) and electrons (e) and assume elastic collisions and no chemical 
reactions in the system. A more complex composition of the plasma, to include other atoms, 
has no relevance for our study, since one can always define an average neutral or ionized 
atom. The equations of conservation of mass, momentum and energy of the three species 
can be written as follows: 
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where all t he notations are standard. Similar system of equation is used in e.g. , iGoedbloed &: Poedts 
( 120041 ) and IZaqarashvili et al.l (120111 ) . The terms R a (a = e,i,n) for elastic collisions can be 

expressed as: 



Re = -p e [Vei(u e - Ui) + V en (u e - U n )} 
Ri = ~ Pi[Vie{Ui ~ W e ) + U in (Ui - U n )} 
Rn = -Pn[Vni{u n ~ Ui) + V ne (u n - U e )] 
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the s um of all three being zero. Here u a p are collisional frequencies. We take these frequencies 
from Spitzerl (1962 ). for collisions between the neutrals and charged particles, and from 
Braginskii Jl965h . for collisions between the charged particles: 



n P e 4 A 



8k B T 

8k^T 
irm er 



-(J e . 



where m in = mim n / (m 
are Oi n = 5 x 10 -19 m 2 and a, 



v„. = a,, \ i <r,„ (5) 

(6) 

m n ) and m en = m e m n / (m e + m n ). The respective cross sections 
10~ 19 m 2 . A is the Coulomb logarithm: 



3/2 



A = 23.4 - 1.15 log 10 n e + 3.45 log 10 T (8) 
with n e expressed in cgs units and T in eV. 

We further assume that the pressure te nsor p can be approx imated by the scalar pressure 
and that the heat flux q is zero. Following iBittencourtl ( 119861 ) . we added up the equations 
for the three species. Under no further approximations, this leads to a system of quasi-MHD 
equations of the form: 
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Fig. 1. — Diffusion coefficients from Eq. (|17[) calculated in the model atmosphere representing a 2nd-order 



thin magnetic flux tube (jKhomenko et al.ll2008at iPneuman et al.lll986 ). 
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is the diffusion velocity, i.e., the difference between the center of mass velocity of the system, 
u, and the velocity of the individual components, u a . We further consider that the diffusion 
velocities are small and neglect the term containing w\ in the definition of the total scalar 
pressure. 



To close the system, the generalized Ohm's law is n eeded to describe the evolution of the 
currents. We use here the form of the Ohm's law from iBraginskiil (119651 ) for slow processes: 
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Fig. 2. — Temperature increment with height estimated after Eq. (|25p in the flux tube model atmosphere, 
assuming only 1% of the existing magnetic field is dissipated (red line). Black line shows the temperature 
run of the VAL-C model atmosphere. 



E + UX B 



Jx B 



(17) 



where J± is the component of current perpendicular to the magnetic field. This form of the 
generalized Ohm's law neglects the temporal variations of the relative ion-neutral velocity, 
the effects on the currents by the partial pressure gradients of the three species, and the 
gravity force acting on electrons. The diffusion coefficients are given by the formulae: 
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(18) 
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(20) 



The three terms on the right-hand side of the generalized Ohm's law are, in order, the classi- 
cal ohmic term, the Hall term and the ambipolar or neutral term. The diffusion coefficients 
defined above depend on the collisional frequency, on the magnetic field strength and on the 
neutral fraction. It is important to note that, for a fully ionized plasma, the parameters p n , 
Vi n and v en are zero. The squared exponent of p n makes tja vanish in this case, due to the 
linear dependence of the collisional frequencies with the density of neutrals (see Eqs. [5] and 
[H]). To have an idea of the expected range of values of these coefficien ts, we have calculated 



them in a model atmosphere representing a 2nd-order thin flux tube (jPneuman et al.lll986 
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Khomenko et al.ll2008al ). The field strength in this model decreases with height, from about 
750 G in the photosp here to 37 G in the chromosphere, and the temperature structure is 
the same as VAL-C (IVernazza et al.lll98ll ). The results are given in Fig. [TJ The classi- 
cal ohmic term reaches its largest values at photospheric heights, between and 500 km. 
However, even there, the Hall term is about one order of magnitude larger than the ohmic 
term. The ambipolar term becomes dominant over the other two from 900 km upwards. At 
chromospheric heights, this term is up to 5 orders of magnitude larger than the ohmic term. 



Using the Ohm's law, the equation for internal energy takes the shape: 



Dp 
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pVu = r]fi J 2 + VAUoJj 



(21) 



As can be seen, the Hall term of the Ohm's law does not appear in the energy equation 
and, consequently, has no impact on the thermal evolution of the system. This is no surprise, 
given that tjh does not depend on the collisions among the species (see Eq. [19]). Since 
we are interested in mechanisms that may lead to the heating of magnetic structures at 
chromospheric heights, we will concentrate in the rest of the paper on the action of the first 
and third terms of the Ohm's law and the Hall term will be removed from Eq. [TTJ 

The induction equation for the temporal variations of the magnetic field is derived from 
the Ohm's law and Maxwell equations. It has the following form: 
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Finally, combining the equations of momentum (Eq. [TO]) , internal energy (Eq. |2~T]) and the 
induction equation (Eq. [22]) . we obtain the following equation for the variations of the total 
energy: 



de tot -> (f B 2 



B(uB) 



(23) 



pug + V B x (77 + r] A )J 



where e tot = pu 2 /2 + £ 2 /2/j + P /(l - 1). 

Here, it is in order to note that the action of the ambipolar (neutral) term in the energy 
equation will stop in two cases. Firstly, when the plasma becomes fully ionized, since tja = 0. 
And, secondly, when the magnetic field evolves to a force-free configuration (J± = 0), since 
J x B = 0. 
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3. Order of magnitude estimates 

The Joule dissipation of currents allows to transform the energy stored in the magnetic 
field into thermal energy. Assuming a scenario where the kinetic energy variations are 
negligible, we can make an order of magnitude estimate of the amount of heating that can 
be achieved by dissipating a given amount of magnetic energy and fully converting it into 
thermal energy. One can write that the loss of magnetic energy is equal to the increase of 
thermal energy, which, to first order approximation, can be expressed as 

^ = At- (24) 

If we further assume that the perturbation of the density is negligible and the totality of 
the pressure perturbations is caused by the temperature variation (leading to the maximum 
temperature increase by magnetic energy dissipation), we can write Ap = pRAT and 

h - 1) 

AT = ^ J -B AB (25) 

pRp 

This approximate equation gives us the amount of temperature increase, AT, when the 
magnetic field is decreased by an amount AB. As an example of how large this energy 
conversion can be, Fig. Ogives the amount of the temperature increase AT according to Eq. 
( |25l) . assuming that as little as 1% of the total magnetic field at a given height is dissipated. 
In this calculation, we used the same thermal and magnetic parameters of the flux tube 
model as in Fig. [TJ We can compare the temperature increase AT (red curve) with the 
actual temperature at each height given by the VAL-C atmospheric model (black curve). It 
shows that an important temperature increase can be reached at heights above 1200-1300 
km, already in the lower chromosphere. Although the assumed 1% is totally unmotivated, 
this example shows that a very small conversion of magnetic energy, even below detectable 
limits in terms of field strength, can lead to a significant amount of heating. 

The time scales at which this dissipation happens depend on the value of the ambipolar 
diffusion coefficient and on the typical scales of the system. Assuming typical velocities of 
motion of V = 10 4 m s" 1 , spatial scales of L = 10 5 m and r] A = 10 7 - 10 8 (at 1500-2000 km, 
see Fig. [TJ the order of magnitude estimate of the time scale gives: 

t « L 2 /r] A = 100 - 1000 sec. (26) 

This time scale is much shorter than in the case when only ohmic diffusion is considered 
(t ~ 10 7 sec, i.e about 4 months). Thus, we may expect that important heating can be 
reached due to ambipolar diffusion in a relatively short time interval of the order of minutes. 
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Fig. 3. — Variations of the magnetic field (black dashed line) and temperature (red solid line) with respect 
to their equilibrium values across the horizontal cut through the vertical magnetic flux tube 2000 sec after 
the start of the simulation. 



The same applies to the Reynolds number. The Reynolds number defined by the am- 
bipolar diffusion is estimated as: 

R m « VL/ VA = 10 - 100, (27) 
while the Reynolds number for the ohmic diffusion in the same conditions is much larger, 



4. Unstratified atmosphere 



In this section, we consider the simple problem of an isothermal atmosphere without 
gravity. We assume the vertical magnetic field depending only on horizontal coordinate as: 



B z (x) = B exp 



2a 2 



(2f 



Gas pressures are prescribed to keep the model in magnetostatic equilibrium. The parameters 
(r = 40 km and B is defined later. 

We solve numerically the equations of conservation of mass (Eq.|UJ), momentum (Eq. [TO]) , 
total energy (Eq. 155]) . and the induction equation (Eq. 155]) . after subtracting the equilib- 
rium conditions from them. These equati on are solved by means of our code mancha 
( iKhomenko et al.ll2008bt Felipe et al.l I201CH ) . The code, as described in the above papers, 
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Fig. 4. — Time evolution of the variations of the internal energy (A^int = -Bint — -E'int(O), black solid 
line); magnetic energy (A_Emag = -Emag — -Emag(O), blue dashed line) and kinetic energy (AE'kin = 
£«in — -Ekin(O), red dotted line), in the simple unstratified model atmosphere containing a vertical flux 
tube. The energies are averaged over the horizontal direction. The variations are shown in per cents from 
the total (horizontally averaged) initial energy of the system £tot(0) = Bint(0) + -Emag(O) + -Ekin(O). 

solves the ideal MHD equations for non-linear perturbations to the magneto-static equilib- 
rium, without physical diffusive terms. We modified the code to include the ohmic and 
ambipolar diffusion terms in the equation of energy conservation and in the induction equa- 
tion, as described above. As our code solves the equations for non-linear perturbations, we 
treat the diffusion terms as perturbations. Without introducing this (or any other) perturba- 
tion no evolution is possible. The equations are solved in two spatial dimensions, though the 
vector quantities are allowed to have three dimensions (2.5D approximation). As the tem- 
perature, and the ionization state, varies with time, we recalculate the ionization balance 
of the atmosphere at each time step, assuming Local Thermodynamic Equilibrium (Saha 
equations). We then update the neutral fraction, p n /p, needed for the calculation of the 
ambipolar diffusion coefficient (Eq. [20]) . 

As initial atmospheric parameters far from the axis of the magnetic structure defined 
by Eq. [28] , we take values of the temperature, pressure and density from VAL-C model at 
1600 km. The value of the magnetic field at the center of the flux tube is limited by the 
condition of the horizontal pressure equilibrium (P mag + P gas = const). According to this 
condition, we set B = 18 G. 

The simple atmosphere considered in this section is initially in equilibrium, obtained 
without considering the diffusion terms. Without external perturbation, it does not evolve. 
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Fig. 5. — Temperature increase AT reached after the 2000 sec of the simulation time in the simple 
unstratified model atmosphere containing a vertical flux tube (red asterisk). The different points are for 
the four different simulations with thermal parameters taken from the VAL-C model atmosphere at the 
corresponding heights. The initial field strength is taken to be Bq = 1500, 100, 18 and 3 G at 0, 600, 
1000 and 1600 km, respectively. This field has been decreased by 9, 5, 2 and 0.2 G after the 2000 of the 
simulations. The black diamonds show the values of AT obtained from Eq. [25j 

But, after introducing the perturbation in the form of diffusion terms, we perturb the initial 
magnetic field structure via the induction equation. Then, as time evolves, this perturbation 
translates to the rest of the variables of the system and it starts to change. Figure |3] shows 
the evolved system 2000 seconds after the introduction of the perturbation. At this time 
moment, the temperature at the center of the flux tube has increased by about 2000 K and 
the magnetic field has decreased by about 2 G. 

The time variation of the different contributions to the total energy are given in Fig. HJ 
It gives the internal, magnetic and kinetic energies, averaged over the horizontal direction 
of the simulation domain. The average kinetic energy is negligible compared to the other 
two and no important plasma motions are produced. The internal energy increases with 
time, and the magnetic energy decreases, in exactly the same amount. Thus, in this simple 
experiment we convert all the released magnetic energy into internal energy by means of 
Joule dissipation of electric currents enhanced due to neutral diffusion (the contribution of 
the ohmic diffusion is very small). As no energy dissipation mechanisms are considered in 
this toy simulation, the internal energy (and temperature) would be expected to ever increase 
with time. As was noted at the end of Section 2, the action of the heating term would stop 
in two cases. In one possible scenario, once sufficiently high temperatures are reached, the 
plasma becomes totally ionized and the ambipolar diffusion coefficient becomes t]a = 0. In 
the second scenario, the magnetic field evolves to a force-free configuration J± = and 
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dissipates a small or large fraction of its energy, depending on the magnetic configuration. 
None of these two conditions are yet reached after 2000 sec of the simulation time shown in 
Fig. [Hand the internal energy shows a tendency to increase with time. 

Obviously, this calculation is just an illustration to show the efficiency of the proposed 
mechanism for chromospheric heating. In a more realistic situation some energy dissipation 
mechanism would balance the heating so that the temperature would not increase unbounded 
to ionize the whole chromosphere. The most probable mechanism to balance the heating 
comes from radiative energ y losses , show n to be a key dissipation mechanism for the chro- 



mosphere since the work of lAthayl (119661 ) 



We have repeated the above experiment for different thermodynamic and magnetic 
field parameters taken at different heights of the VAL-C model atmosphere. The resulting 
temperature increase reached after 2000 sec of the simulation time is shown in Fig. For 
comparison, we show the values of AT obtained from the simple formula (Eq. 125]) . The most 
important heating is reached at 1000 km height. The values given by the formula are in 
order-of-magnitude agreement with those obtained in the numerical experiment, suggesting 
that this formula can be used as an estimation. As the amount of heating depends on the 
initial magnetic field, less heating is reached at 1600 km, where we took the initial magnetic 
field to be only B = 3 G. 



5. Vertical flux tubes 

In this section, we describe the results of a similar simulation, but now in a gravitation- 
ally stratified isothermal model atmosphere. As initial magnetic field structure, we take the 
vertical flux tube with B z (x) given by Eq. ff28|) . We made five simulations with the value 
of Bq varying from 300 to 900 G. The atmosphere has initially a constant temperature of 
T = 5179 K and both pressure and density are stratified with a scale height of 157 km. Note 
that now the atmosphere is not in horizontal force balance since B z is constant with height. 
The size of the simulation box is 1 Mm by 1.65 Mm with a resolution of 10 km. We take 
periodic boundary conditions in horizontal direction and closed boundary in both vertical 
directions. Horizontal motions of the magnetic field lines at the top and bottom boundary 
are allowed. 

As can be appreciated from Fig. [TJ the ambipolar diffusion coefficient t]a increases 
exponentially with height. The diffusion term in the energy equation is proportional to rjA- 
Thus, we expect the largest heating at the upper part of the simulation domain. To avoid 
any possible numerical artifact due to the conditions at the upper boundary we artificially 
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decreased t)a gradually to zero in the upper 15 grid points (150 km) of the domain. To that 
aim, we used a third-order polynomial that matches the value of tja and its derivative in one 
end and reduces down to zero, and zero derivative, at a distance of 7 grid points (70 km). 
The rest of the points, up to 150 km, have tja equal to zero. This way, we make sure that 
the heating starts inside the physical domain, not influenced by the upper boundary. 

Figure [6] shows the time evolution of the temperature and magnetic field strength in the 
simulation with initial B Q = 300 G. The heating starts immediately at the central part of 
the tube in a few seconds time scale. The appearance of the temperature perturbation at 
the upper part of the simulation domain goes in agreement with the current distribution in 
the initial magnetic field configuration and also with the fact that the ambipolar diffusion 
coefficient reaches maximum values at 1.5 Mm according to our simulation setup, and then 
reduced to zero at the topmost layers. It is evident from the figure that, in about 10-seconds 
time, the temperature increases about 2000 K from its initial value of To = 5179 K at heights 
around 1.5 Mm, with a much lower variation in those layers where t\a is artificially decreased. 
It must be noted that the perturbation introduced by tja also causes plasma motions. At the 
initial stage of the simulation, these motions are upflows at the central part of the tube. Due 
to the action of the advective term W(ue to t), the velocity motions transport hot material to 
the upper layers and heat them. This heating is thus a consequence of the transport of hot 
plasma and not the direct consequence of ambipolar diffusion. The amplitudes of the plasma 
motion do not exceed 100 m/sec and decrease with time as the tube evolves. The diffusion of 
the field causes its "opening" with height. The field becomes progressively more horizontal 
with time at heights above 0.8 Mm, which is also a direct consequence of the increase of t)a 
with height. The expansion of the field, together with the outward plasma motions, cause 
the temperature increase also in the surroundings of the flux tube, since the temperature 
perturbation is transported by the advective term of the energy equation. But the heating is 
the largest in the central part of the tube where the magnetic field is the strongest. Also the 
time scales of heating inside and outside the flux tube are very different. While the central 
part of the tube is heated by several kK in a time interval of a few seconds, the surroundings 
take longer time to evolve and are heated by several hundred K in 5-10 minutes. However, 
the details of the behavior of the flux tube surroundings are not relevant for our study. 

Another visual impression from Fig. [H] is that the heating at the central part of the 
tube stops after 100 sec of the evolution. It actually stops sooner, but this time moment is 
not shown in the figure. The explanation for that is provided in Fig. [71 This figure shows 
the time evolution of the ionization fraction, the ambipolar diffusion coefficient T]a and the 
characteristic time scale given by the ambipolar diffusion t ~ L 2 /t]a at the central part of the 
flux tube at a height of 1.4 Mm. During the first 20 seconds of the evolution, the changes are 
very significant. The material gets heated and the ionization fraction increases from 10 -3 to 
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0.5. This produces an immediate decrease of the ambipolar diffusion coefficient by about 3 
orders of magnitude. Thus, the time scale also increases abruptly. Initially the heating takes 
place on about few tens of seconds time scale, but after first 40 seconds the characteristic 
time becomes about 20000 sec (about 5 hours). This explains the slowness of the evolution 
of the temperature at the central part of the tube after the first few tens of seconds. 

Finally, Figure [8] gives the temperature variation as a function of time at the center of 
the flux tube, averaged between 1 and 1.5 Mm heights, for all five simulations with different 
Bq. As expected, since t/a depends on the magnetic field, the temperature increase is more 
significant for the larger fields. AT reaches a maximum of 3.5 kG for Bq = 900 G at t = 70 
sec, making T = 8500 K. The magnetic field strength has been decreased by only about 3 G. 
In all cases, the temperature rise is extremely quick during the first few tens of seconds. It 
then reaches a stationary stage after t = 100 — 200 seconds, depending on the field strength. 
In the stationary stage, the values of the temperature stay between 7000 and 8000 K. These 
values are larger tha n those given by the VAL-C model atmosphere at 1.5 Mm (6400 K, see 



Vernazza et al.lll98ll ). 



The above simulations have several concerns. Firstly, the vertical flux tube is initially 
not in horizontal force equilibrium, which might cause an excess of the horizontal plasma 
motions and other side effects. In this respect, we would like to remind that our code does not 
solve the full differential equations. After subtracting the equilibrium conditions from Eqs. HI 
[TUl I2"3"l and [221 the resulting differential equations for the evolution of the perturbations are 
solved. This means that, even if our flux tube is not initially in equilibrium, it cannot 
evolve if no perturbation is applied to it. In our case, the ambipolar diffusion represents this 
perturbation, and without it no evolution is present. Thus, the results of our calculation 
are determined by this parameter and not by the deviations from equilibrium. We can not 
disregard, though, that the exact evolution (especially of the velocities) may be influenced 
by the non-perturbed values. These cannot, however, have any influence on the thermal 
evolution, which is fully determined by the ambipolar diffusion. 

The second concern may come from the adopted magnetic field strength, which is, 
possibly, too large as for a quiet chromosphe ric region. Measuring magnetic fields in the quiet 



chromosphere may be challenging (see, e.g.. lManso Sainz &: Trujillo Buenoll2010l ) and actual 



measurements are lacking. In active regions, where the polarization signal s are stronger, 



non-LTE in versions of the infrared Ca II triplet and Fe I lines have allowed ISocas-Navarro 



(120051 . 120071 ) to constrain the magnetic field strength to about 1.5 kG above sunspots and 
about 1 kG in a network region. As far as we are aware, the measurements for the quiet 
chromosphere are still missing. But, as a reasonable choice, one may assume values of about 
tens of Gauss. 



- 16 - 



To address these two concerns, in the next section we describe a simulation done in a 
more realistic flux tube model in equilibrium, with a lower field strength at chromospheric 
level, more representative a quiet area. 



Second-order thin flux tube 



To simulate non-active chromosp heric conditions, we used the 2nd-order thin m agnetic 
flux tube as initial model atmosphere (jPneuman et al. 1986 ; Khomenko et al. 2008a ). Previ- 



ously, the same flux tube model was used for simulations of wave propagatio n from the photo 



sphere to the chromosphere and it is described in the corresponding papers (IKhomenko et al. 



2008al ; Khomenko et al.ll2008bl ). The model represents a series of flux tubes that merge af- 
ter some height in the chromosphere, preventing them from excessive opening with height. 
Thus, the magnetic field configuration is non-potential and there are currents in the sys- 
tem. Figure IH] illustrates the initial distribution of some parameters of the flux tube model. 
The magnetic field strength (Fig. [HJ upper left) has both vertical and horizontal gradients. 
The colors of the figure are saturated at the bottom part of the structure to highlight the 
variations of the magnetic field in the chromosphere. The field is stronger at the central 
part of the tube, with about 37 G from 1 Mm upwards. The ionization fraction defined 
as p e / p (Fig. EH upper right) varies in agreement with the temperature distribution of the 
VAL-C model atmosphere. It has its lowest values in the photosphere, dropping down to 
10~ 4 and reaching about 0.7 in the upper part of the domain, at 1.5 Mm. The ambipolar 
diffusion coefficient (Fig. [9j bottom left) changes from 10 2 in the photosphere up to 10 s in 
the chromosphere inside the flux tube (note its artificial decrease in the upper 15 grid points, 
120 km, done for the same reasons as in the simulation shown in the previous section). Out- 
side the flux tube, its value is very small, as the magnetic field strength outside the tube is 
negligible (we set magnetic field outside the flux tube equal to some small value for reasons 
of numerical stability of the simulations). The most interesting quantity is the one shown 
at the bottom right panel of Fig. ED It gives a proxy of the diffusion term calculated as the 
square root of the current J 2 multiplied by the ambipolar diffusion coefficient (see Eq. |2"T|) . 
As the variations of the magnetic field are stronger near the borders of each individual tube, 
the current is also stronger there. The absolute value of J 2 decreases according to the drop of 
the magnetic field strength with height, while tja exponentially increases. When multiplied 
one by the other, the quantity J 2 r\A has largest values at the upper part of the domain above 
1 — 1.3 Mm, near the borders of the tubes. This term is responsible for the heating, therefore 
we expect to observe the largest heating of the flux tube atmosphere at the locations where 
J 2 i]a is important. 
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Figure [TO] gives several snapshots of the evolution of temperature in the flux tube after 
the diffusion terms were introduced as perturbation. Same as before, the boundary condi- 
tions are taken periodic in the horizontal direction (equivalent to an infinite series of flux 
tubes); the domain is closed at the top and bottom boundaries; and t/a is decreased smoothly 
to zero at the upper 100 km of the simulation domain. Except for velocity variations (unim- 
portant for the conclusions of the present study), closed boundary conditions do not affect 
the variations of any other quantity in our simulations. 

Similar to the case of vertical flux tubes (described in Sect. EJ), the heating starts at the 
upper part of the domain, close to the tube borders. This behavior is expected because the 
term responsible for the heating (~ J 2 i]a, see Fig. [9]) is orders of magnitude larger at these 
locations. The upper 120 km of the atmosphere, with artificially reduced t]a, are not heated 
during the first few seconds, although the heat is transported there by the advective term 
in the energy equation from the material just heated below. The heating gradually extends 
beyond the tube borders and extends down into the atmosphere. The temperature at 1900 
km reaches 8800 K after 800 sec of the simulation, which is about 2000 K above its initial 
value. The time scales of the heating are longer compared to the case of the vertical tubes 
with larger field strength. There is still no saturation after 800 sec of the simulation and the 
temperature is increasing slowly with time. 

The velocity field developed in this simulation is different from the case of vertical flux 
tubes (Sect. E}. The motion is mostly vertical and periodic in time, changing from upflows 
to downflows. Such motion is a direct consequence of the closed upper and lower boundary 
conditions and the periodicity of our magnetic structure (series of flux tubes in horizontal 
direction with mostly vertical field in the chromosphere). There is also no expansion of the 
field with height. The expansion is not possible as the tubes are located one next to the 
other. 

Temporal variations of the ionization fraction, ambipolar diffusion coefficient and the 
parameter L 2 /t)a are shown in Fig. [TTJ The values are taken at 1900 km close to the border 
of the tube at horizontal location X = 0.6 Mm. The initial ionization fraction at this 
height is already rather high. During the first 50 seconds of the simulation, it increases from 
0.7 to about 0.9. Accordingly, the ambipolar diffusion coefficient drops by about an order 
of magnitude after the first 50 seconds. Later, the evolution slows down and the further 
decrease of t]a is more gradual. The change of the characteristic time scale is shown in the 
right panel of Fig. [TJJ The time scale increases from 50 seconds (initially) up to 300-400 
seconds (after 400 sec of the simulation). 

In the case of the vertical tubes (Sect. [5]), the variations of t ~ L 2 /i]a were stronger 
and the values of t at the later stages of the evolution were much larger. The difference 
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appears because, firstly, the initial magnetic field strength of the vertical tubes was much 
larger, causing a rapid heating and ionization within the first 20 seconds of the simulation. 
Secondly, we considered lower heights up to 1500 km. The value of tja changes more than 
one order of magnitude between 1500 and 2000 km (Fig. [TJ and so does the time scale. 
Thus, while the evolution gets slower in time in the simulation with a 2nd order flux tube, 
it is not as slow as in the case of the vertical tubes. After the initial phase of the rapid 
evolution, the temperature in our 2nd order flux tube continues increasing on scales about 
5-6 minutes. Obviously, in a more realistic modeling, this constant energy input is expected 
to be balanced by some dissipation mechanism. Otherwise, the simulation would eventually 
become unstable. However, here we are interested in understanding and evaluating the action 
of the heating (not cooling) term. For that reason, we considered a relatively simple situation 
that allows us to better constrain the physics. Clearly, the final judgment about the heating 
efficiency of the proposed mechanism will require the i nclusion of ra diative energy exchange, 
a key mechanism for the chromospheric energy losses (lAthayl 119661 ). 

Finally, Figure [12] gives the temperature as a function of height at the horizontal location 
X = 0.6 Mm, close to the tube border, at different time moments. The atmosphere gets 
heated above 1100 km and the amount of temperature increase is larger at larger heights, 
except in the uppermost layers at the beginning of the evolution, where the ambipolar 
diffusion does not operate. The magnetic field strength has decreased on average by about 
0.01 G at heights 1100-1900 km and X = 0.6. Thus, we have dissipated only about 0.03% 
of the existing magnetic field. Fig. [T2] can be compared with Fig. [2] where we calculated the 
temperature increase AT assuming the same model atmosphere, but with 1% of the existing 
magnetic field at each height being dissipated. Given that we have dissipated 30 time less 
magnetic field than it was assumed for the calculations of AT in Fig. |2j the results of our 
numerical simulations are in order of magnitude agreement with the simple estimation by 
Eq.| 



7. Conclusions 

In this work, we have investigated the heating of the solar magnetized chromospheric 
plasma produced by current dissipation (Joule heating). The current dissipation is enhanced 
by orders of magnitude due to the presence of a significant amount of neutral atoms in the 
chromospheric plasma not entirely coupled by collisions. The net relative motion between 
the ionized and neutral plasma components results in an additional diffusion term (ambipolar 
or neutral diffusion). The value of the ambipolar diffusion coefficient depends on the fraction 
of neutral atoms and on the magnetic field strength and can be as much as five orders of 
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magnitude larger in the chromosphere than the classical ohmic diffusion coefficient. This 
implies evolution time scales of about 10-100 seconds. Thus, all the processes related to the 
ambipolar diffusion, including the plasma heating, happen very quickly. 

We have performed analytical estimates and numerical simulations in different models of 
the amount of plasma heating due to diffusion by neutral particles. The simulations covered 
a variety of scenarios, from a simple atmosphere without gravity to a more complex flux 
tube model, resembling a chromospheric quiet network area. All our calculations show that 
the lower chromosphere can be heated on temporal scales of 10 2 seconds, giving rise to a 
temperature increase of several kK. The only requirements for the heating is the existence of 
a non force-free magnetic field. We obtain that the temperature increase is significant above 
a height of 1 Mm, depending on the strength of the magnetic field. 

From those results, we conclude that current dissipation enhanced by the action of 
ambipolar diffusion is an important process that is able to provide a significant energy input 
into the chromosphere. We suggest that this mechanism should be included into the future 
models of chromospheric heating. Here we would like to stress that, apart from possible 
dissipation mechanisms, the action of the ambipolar term in the energy equation would stop 
in two cases. Firstly, this will take place if the temperature increase is so significant that all 
atoms become ionized leading to t\a = 0. Obviously, this is not the case of the chromosphere 
which is observed to be, at least, partially neutral. Note that our calculations indicate that 
the time scales of heating may become very long after a fast initial temperature increase (see 
Figures 171 and [TTT) . This would mean that the actual complete ionization is never reached 
in a reasonable amount of time. Secondly, the heating due to ambipolar diffusion will stop 
if the magnetic field configuration becomes force-free, meaning J \\ B, i.e. J± = 0. This 
might lead to different temperatures and ionization fractions for different starting magnetic 
field configurations and will need further detailed investigation. Finally, the action of the 
ambipolar diffusion term has to be balanced by a dissipation mech anism. It is well known 



since the 60's that the chromosphere is cooled by radiative losses (lAthaylll966l ). Whether 
the proposed mechanism is able to balance these losses, and whether it is able to provide 
sufficient energy to support the chromospheric high temperature, has to be investigated in 
future realistic 3D MHD chromospheric simulations, including a detailed radiative transfer 
and all complex chromospheric physics. 

Here, we would just like to give hints that the temporal scales of the heating are defined 
by the magnetic field strength via the ambipolar diffusion coefficient t]a (Eq. [20]) . For the 
stronger fields considered by us (B = 300 — 900 G), the heating happens very quickly, 
achieving a temperature increase of 2-3.5 kK in a 1 minute time. For the weaker field 
strengths (30-40 G), expected in the quiet chromosphere, the plasma is heated by 2 kK in 
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a 10 minutes time. These time scales are comparable to the scales of radiative cooling of 
the plasma (a simple Newton fo rmula gives characteristic cooling times at this height of the 
order of 10 2 — 10 3 seconds, see iMihalas &: Toomrdll982l . figure 3). This argument may be 
used to suggest that one process might, in fact, balance the other. 

The radiative damping would tend to decrease the temperature in the chromosphere. 
This, in turn, would produce a decrease of the ionization fraction, a subsequent increase of tja 
and a decrease of the heating time scale (see Figs. [7] and [TT]). The neutral diffusion term would 
be put in action again, leading to a new temperature enhancement. The detailed temperature 
balance would depend on the details of non-LTE chromospheric radiative transfer and on 
the scales of the non-LTE ionization equilibrium. 

The present study has several limitations. The calculation of the ionization balance is 
done in LTE conditions, assuming instantaneous ionization — recombination. Several recent 
studies have shown that the ionization-recombination of (at least) hydrogen and calcium in 
the solar chromosphere does not follow LTE relations or instantaneous stati stical equilibrium 
(ILeenaarts fc Wedemeyer-Bohmll2006l ; IWedemeyer-Bohm fc Carlssonll201ll ). In addition, we 
have discussed the energy losses by radiative damping, but have not included them in our 
simulations. Thus, in the simulations shown in this paper, the temperature would ever 
increase until the chromosphere is (unrealistically) completely ionized or until the magnetic 
field configuration becomes force-free. This, however, was done on purpose since we wish 
here to investigate the action of a newly proposed mechanism in simple situations, to better 
understand its implications into the general picture. In our future studies, we plan to include 
the effects of radiative damping in our calculations. 
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Fig. 6. — Time evolution of temperature in an (initially) vertical magnetic flux tube described in Sect. [5] 
The color scale is given on the vertical bars, same for all panels. The background grey color corresponds 
to a temperature of 5179 K on the color bar. Vertical white lines are magnetic field lines. Arrows 
show the velocity field. The initial magnetic field strength is Bo = 300 G. The time after the 
start of the simulation is indicated on each panel. Two horizontal blue lines at the top mark 
the boundaries of the regions where tja is decreased to zero (lower line) and where tja — 
(upper line). 
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Fig. 7. — Time variation of the ionization fraction p e /p (left); ambipolar diffusion coefficient t/a (middle); 
and characteristic time scale L 2 /tja (right) at the height of 1500 km at the flux tube center in the simulation 
with B Q = 300 G from Sect. 
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Fig. 8. — Time variation of the average temperature at heights 1000-1500 km at the center of the flux 
tube in the simulations with B = 300, 450, 600, 750 and 900 G from Sect. E3 
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Fig. 9. — Initial distribution of the magnetic field (upper left); ionization fraction p e /p (upper right); 
ambipolar diffusion coefficient tja (bottom left); and the quantity J 2 rjA (bottom right) in the 2nd order flux 
tube model. Note that t/a is decreased artificially to zero at the upper 120 km (15 grid points). 
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Fig. 10. — Same as Fig. [6] but for the 2nd order thin flux tube. For better visual comparison we keep 
unchanged the temperature structure at the left hand side of the domain (from to 0.42 Mm), though the 
variations are present in the simulations. 




Fig. 11. — Same as Fig. [7] but for the 2nd order flux tube model at horizontal distance 0.6 Mm and height 
1900 km. 
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Fig. 12. — Temperature as a function of height at horizontal distance 0.6 Mm of the 2nd order flux tube 
from Fig. 1101 Different lines are separated 20 sec in time with progressively more red colors indicating larger 
times till 800 sec since the start of the simulation. Note that the evolution gets slower with time. Two 
horizontal blue lines mark the boundaries of the regions where tja is decreased to zero (1880 km) and where 
riA = (1940 km). 



